Molecular Simulations and Drug Discovery of Adenosine Receptors

G protein-coupled receptors (GPCRs) represent the largest family of human membrane proteins. Four subtypes of adenosine receptors (ARs), the A1AR, A2AAR, A2BAR and A3AR, each with a unique pharmacological profile and distribution within the tissues in the human body, mediate many physiological functions and serve as critical drug targets for treating numerous human diseases including cancer, neuropathic pain, cardiac ischemia, stroke and diabetes. The A1AR and A3AR preferentially couple to the Gi/o proteins, while the A2AAR and A2BAR prefer coupling to the Gs proteins. Adenosine receptors were the first subclass of GPCRs that had experimental structures determined in complex with distinct G proteins. Here, we will review recent studies in molecular simulations and computer-aided drug discovery of the adenosine receptors and also highlight their future research opportunities.


Introduction
Adenosine (ADO) is an endogenous nucleoside, which regulates multiple biological functions by activating adenosine receptors (ARs) [1][2][3]. ARs belong to the class A G protein-coupled receptors (GPCRs), representing primary targets of approximately 1/3 of marketed drugs [4]. Four subtypes of ARs are expressed in human bodies, including A 1 AR, A 2A AR, A 2B AR and A 3 AR. The A 1 AR and A 2A AR are high-affinity receptors for ADO, whereas A 2B AR and A 3 AR are low-affinity receptors (Table 1) [5,6].
During function, the A 1 AR preferentially couples to the G i/o proteins to inhibit the production of cAMP by regulating the activity of adenylyl cyclase (AC) [3,7]. It could regulate many cellular responses, such as inhibition of Ca 2+ conductance and stimulation of phospholipase C, K + conductance, phosphoinositide 3 kinase (PI3K) and mitogen-activated protein kinase (MAPK) [8,9] (Table 1). The A 1 AR is thus considered as a potential drug target for treating myocardial ischemia [10], cardiovascular disorders [11,12], obesity [13] and cancers [14,15]. Particularly, one agonist and four antagonists of A 1 AR were approved in the market ( Table 1).
The A 2A AR prefers to couple to the G s protein to activate AC, resulting in the activation of cAMP-dependent protein kinase A (PKA), protein kinase C (PKC), MAPK and ion channels (Table 1) [16]. The A 2A AR is widely considered as a potential drug target for treating cardiovascular disorders [11,12], obesity [13], Parkinson's disease (PD) [17,18] and cancers [14,15]. Particularly, a selective A 2A AR antagonist Istradefylline was approved for the treatment of PD [19].
The A 2B AR binds with the G s and G q proteins to induce the PKA signaling to increase the level of cAMP and stimulate phospholipase C and MAPK [20]. The A 2B AR is recog- The A 3 AR couples to the G i protein to inhibit AC and decrease cAMP accumulation and PKA activity. Additionally, the A 3 AR could bind with the G q protein to stimulate phospholipase C, resulting in increased Ca 2+ levels and modulation of PKC activity [25]. In addition, they may activate the phospholipase C pathway through the G βγ subunit [26]. Increasing attention is paid to the A 3 AR for the drug development against inflammation [27], glaucoma [28], rheumatoid arthritis [29] and stroke [30].
Advances of techniques in X-ray crystallography, cryo-electron microscopy (cryo-EM) and Nuclear Magnetic Resonance (NMR) have generated a number of structures for GPCRs including the ARs. They have stimulated numerous structure-based and computer-aided drug design (CADD) developments for ARs [31]. To date (January 2022), Protein Data Bank (PDB) [32] has recorded approximately 64 structures of ARs under different functional states (Table 2). However, there has been no experimental structure for A 2B AR and A 3 AR yet. Additionally, only static snapshots of ARs could be captured in the PDB structures. It is still challenging for experimental techniques to probe flexibility of the ARs, which plays a critical role in regulating their biological functions and designing potent and selective drug molecules. Molecular dynamics (MD) simulations [33] have proven useful to model protein dynamics and ligand binding/dissociation processes, which has greatly facilitated the rational drug design targeting ARs. Here, we present a review of simulation studies that have revealed important insights into the dynamic and functional mechanisms of ARs, as well as drug discovery efforts targeting the ARs.

Activation of Adenosine Receptors
The classical view of GPCR activation involves conformational states between the active and inactive states [77]. Agonist binding shifts the conformational ensemble of receptor to the activated state to enable coupling of the G protein. Recent Figure 1). Comparison between the active and inactive structures of the A 1 AR [35,38] revealed a large outward movement of transmembrane helix 6 (TM6) intracellular domain during activation of the receptor, which opened the intracellular pocket for G protein binding (Figure 1). Meanwhile, the receptor activation was accompanied by adjustments of the TM7, helix 8 (H8), extracellular loops (ECLs) and ligand binding pocket (Figure 1). A number of important structural motifs, named "microswitches", play a critical role in the activation of ARs, including the R 3.50 -E 6.30 salt bridge, the W 6.48 "rotamer toggle switch" and the Y 7.53 "tyrosine toggle switch" [78]. The Ballesteros-Weinstein numbering [79] is used for residues in the GPCRs. Although NMR has been used to probe populations of different activation states of the A 2A AR [58,80,81], it remains challenging to probe dynamic conformational transitions among these different states using experimental techniques. In this regard, MD simulations have been widely applied to study the activation mechanism of ARs.
The R 3.50 -E 6.30 salt bridge was very stable in MD simulations of the apo A 2A AR in the inactive state [82,83]. Rotameric transition of the W 6.48 was observed in MD simulations of the agonist 5 -N-carboxamidoadenosine (NECA)-bound A 2A AR, which allowed for water movement through the TM bundle [84]. Such movement further induced a rotational switch of the residue Y 7.53 , inward movement of TM7 on the intracellular side and outward movement of TM6, leading to the opening of the intracellular pocket to facilitate G protein binding [85]. Advances in computing hardware (e.g., GPU and Anton) and software developments have enabled longer conventional MD (cMD) simulations [86]. Even so, cMD is often limited to typically hundreds of nanoseconds to tens of microseconds [87,88]. However, activation of GPCRs including ARs is still beyond the accessible time scale of cMD, which often occurs over milliseconds or even longer timescales [89].
Enhanced sampling techniques that could observe much longer timescale events within shorter simulation time have thus been applied to explore AR activation and deactivation [90,91]. Adaptive sampling combined with Markov State Models (MSMs) allowed for calculation of the activation energy landscape of the apo A 2A AR. The apo A 2A AR sampled four predominant states, including the inactive antagonist bound-like state, the inactive apo intermediate state, the agonist-competent state and the active state [92]. Moreover, Li et al. performed extensive Metadynamics simulations to explore the activation and deactivation mechanism of the A 2A AR [93]. The A 2A AR was simulated in the agonist-bound, antagonist-bound and apo forms. The simulations sampled distinct conformational states of the A 2A AR that resembled the receptor active and inactive states, involving marked conformational transitions of the W 6.48 toggle switch. Three distinct regions in the orthostatic binding pocket were identified to be responsible for the ligand binding affinity, selectivity and agonism/antagonism, respectively. In addition, deactivation of the A 1 AR was observed in our recent Gaussian accelerated molecular dynamics (GaMD) simulations of the active receptor after removing the G protein [39]. In summary, MD simulations, especially enhanced sampling, have provided unprecedented insights into the activation mechanism of ARs at the atomistic level.

Specific G Protein Coupling to Adenosine Receptors
The odd ARs (A 1 AR and A 3 AR) preferentially couple to the G i/o proteins, while the even subtypes (A 2A AR and A 2B AR) preferentially couple to the G s proteins. However, increasing studies suggest that GPCRs including the ARs could couple to multiple G proteins [94][95][96][97]. MD simulations were applied to investigate the dynamic AR-G protein interactions. The role of G protein in stabilizing the active state of A 2A AR was investigated by Lee et al. [98]. Four different conformations of A 2A AR, including the inactive, activeintermediate and fully active in the presence and absence of the mini-G s protein, were used as simulation starting structures. In comparison to the inactive state, lower agonist fluctuations and decreased entropy in the ECLs were observed in the active-intermediate state. In the fully active G protein-bound state, the entropy of ECLs was the highest. The volume of the receptor orthosteric pocket decreased to enable tighter contacts with the  5UEN) and active agonist-G protein-bound (PDB: 6D9H) A1AR. (A) Outward movement of transmembrane helix 6 (TM6) in the active agonist-bound receptor (red) induces opening of the intracellular pocket for binding with G protein as compared to the inactive antagonist-bound conformation of the receptor (blue). (B) The receptor activation is accompanied by adjustments of the ligand binding pocket, extracellular loops (ECLs) and the W 6.48 "rotamer toggle switch". Antagonist DU172 (blue) and agonist adenosine (ADO, red) occupy different regions of the orthosteric pocket in the inactive and active conformations of the A1AR receptor, respectively. (C) "Microswitches" play a critical role in the activation of adenosine receptors, including the R 3.50 -E 6.30 salt bridge and the Y 7.53 "tyrosine toggle switch".
The R 3.50 -E 6.30 salt bridge was very stable in MD simulations of the apo A2AAR in the inactive state [82,83]. Rotameric transition of the W 6.48 was observed in MD simulations of the agonist 5′-N-carboxamidoadenosine (NECA)-bound A2AAR, which allowed for water movement through the TM bundle [84]. Such movement further induced a rotational switch of the residue Y 7.53 , inward movement of TM7 on the intracellular side and outward movement of TM6, leading to the opening of the intracellular pocket to facilitate G protein binding [85]. Advances in computing hardware (e.g., GPU and Anton) and software developments have enabled longer conventional MD (cMD) simulations [86]. Even so, cMD is often limited to typically hundreds of nanoseconds to tens of microseconds [87,88]. The receptor activation is accompanied by adjustments of the ligand binding pocket, extracellular loops (ECLs) and the W 6.48 "rotamer toggle switch". Antagonist DU172 (blue) and agonist adenosine (ADO, red) occupy different regions of the orthosteric pocket in the inactive and active conformations of the A 1 AR receptor, respectively. (C) "Microswitches" play a critical role in the activation of adenosine receptors, including the R 3.50 -E 6.30 salt bridge and the Y 7.53 "tyrosine toggle switch".
In a recent study [99], we employed GaMD simulations on the cryo-EM structures of native agonist ADO-bound A 1 AR-G i and NECA-bound A 2A AR-G s protein complexes [38,52], as well as "decoy" complexes generated by switching the G proteins (A 1 AR-G s and A 2A AR-G i ). GaMD simulations suggested that slight differences of agonist NECA binding in the A 2A AR were observed upon changing the G s protein to the G i ( Figure 2C), while significantly increased fluctuations in the A 1 AR and ADO were identified upon changing the G i protein to the G s ( Figure 2D). The agonist ADO sampled two different binding poses ("L1" and "L2") when the A 1 AR coupled with the G s protein. In the "L2" binding pose, ADO formed interactions with residues Y 1.35 and Y 7.36 in the sub-pocket 2 of the A 1 AR as described earlier [35]. Only one stable low-energy conformation was observed for each of the A 1 AR-G i and A 2A AR-G s complexes as in the cryo-EM structures (Figure 2A,B), being similar for the A 2A AR-G i complex ( Figure 2C). While multiple states of the ADO agonist and G s protein were sampled in the A 1 AR-G s system ( Figure 2D). Simulation results thus indicated that the A 1 AR preferred to couple with the G i protein to the G s ( Figure 2E), while the A 2A AR could couple with both the G s and G i proteins ( Figure 2F), being highly consistent with experimental findings of the ARs [94][95][96]. Further detailed analysis of the simulation trajectories suggested that remarkably complementary residue interactions at the protein interface led to the specific AR-G protein coupling, which involved mainly the receptor TM6 helix, the G α α5 helix and α4-β6 loop. In summary, the GaMD simulations have provided important insights into the dynamic mechanism of specific GPCR-G protein interactions.

Biased Agonism of Adenosine Receptors
Biased agonism is a phenomenon in which binding of different agonists to a target GPCR promotes distinct receptor conformations that bias cellular signaling toward and away from a subset of pathways, e.g., activation of different G protein or arrestin [100][101][102][103]. Biased agonism was first introduced by Jarpe et al. [103] and has been targeted for developing selective therapeutics of GPCRs [102,104]. The first A 1 AR biased agonist LUF5589 was identified by assessing A 1 AR bias with 800 A 1 AR agonists and antagonists in 2013 [105]. The AR biased agonism has been explored by MD simulations. By combining MD simulations, Gαi/o subunitand β-arrestin-specific cellular signaling assays, Wall et al. [106] identified that the A 1 ARselective agonist, BnOCPA, was a biased agonist in exclusively activating the G ob protein. MD simulations were applied on the A 1 AR bound by the biased agonist BnOCPA, neutral agonists ADO and HOCPA and an antagonist PSB36. The simulations suggested that BnOCPA engaged with the same receptor interactions as neutral agonists ADO and HOCPA. However, a distinct partial transition of the N 7.49 PxxY 7.53 backbone from the active to inactive state was observed in one of the BnOCPA-bound A 1 AR simulations. The α5 helix (GαCT) of the G protein (G i2 , G oa , G ob ) was dynamically docked to the HOCPA-and BnOCPA-bound active A 1 AR structures to study the agonist-driven interaction between the A 1 AR and the G protein. MD simulations suggested that the GαCT of G ob docked to the A 1 AR via a metastable state (MS1) relative to the canonical state (CS1). The CS1 corresponded to the canonical arrangement as captured in the cryo-EM structure of the A 1 AR-Gi (PDB: 6D9H), whereas the MS1 was similar to the non-canonical state in the neurotensin receptor, being suggested as an intermediate on the way to the canonical state [107]. In contrast, docking of the GαCT of G oa and G i2 to the A 1 AR formed MS2 and MS3 states. The MS2 was similar to the β 2 -adrenergic receptor-G αs CT complex [108], which was proposed to be an intermediate on the activation pathway and play an important role in G protein specificity. Additional MD simulations were performed on the BnOCPA-and HOCPA-bound A 1 AR in complex with the entire G oa and G ob proteins, respectively. The main differences between the G oa and G ob proteins comprised the formation of transient hydrogen bonds between the α3-β5 and α4-β6 loops of G oa and H8 of the A 1 AR. Overall, G oa interacted more with residues at the TM3 and ICL2 of receptor, while TM5, TM6 and ICL1 were more engaged by G ob . Particularly, residues R 7.56 and I 8.47 showed a different propensity to interact with G oa or G ob proteins. Therefore, a particular A 1 AR conformation stabilized by BnOCPA may favor different intermediate states during the binding process of G oa and G ob proteins.
In previous studies, bitopic ligands that comprise of orthosteric and allosteric ligand binding moieties were shown to exhibit biased agonism in the A 1 AR [100,109]. The supervised MD (SuMD) simulations were performed to capture the binding of bitopic agonist VCP746 (biased agonist) and its allosteric part (VCP171) to the A 1 AR [110]. The VCP746 sampled the most stable binding pose when the orthosteric site was occupied by the adenosine moiety, suggesting that the agonist component of the VCP746 played a particularly important role in the binding. VCP746 interacted with many receptor side chains and the ECL2 backbone atoms. The linker was captured to insert between the E172 ECL2 and M 5.35 side chains in the simulations. The allosteric moiety part VCP171 sampled many orientations on ECL2 but stabilized in a binding mode near ECL2 when the adenosine moiety reached the orthosteric binding site. Side chains E172 ECL2 and K173 ECL2 formed a sort of saddle for the allosteric moiety, which often oriented the 3-(trifluoromethyl)phenyl group toward the hydrophobic pocket formed by residues K173 ECL2 and I167 ECL2 . In comparison, the positive allosteric modulator (PAM) VCP171 usually oriented its 3-(trifluoromethyl)phenyl moiety toward the top of ECL2. Taken together, these computational findings suggested a different binding mode for VCP171 as part of VCP746 and proposed an allosteric site of ECL2 as involved in the observed bias ( Figure 3).

Molecules 2022, 27, x FOR PEER REVIEW 8 of 27
GaMD simulations have provided important insights into the dynamic mechanism of specific GPCR-G protein interactions. Summary of specific AR-G protein interactions: (E) the ADO-bound A1AR prefers to bind the Gi protein to the Gs. The latter could not stabilize agonist ADO binding in the A1AR and tended to dissociate from the receptor. (F) The A2AAR could bind both the Gs and Gi proteins, which adopted distinct conformations in the complexes. Adapted from reference [99] with permission from American Chemical Society. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Biased Agonism of Adenosine Receptors
Biased agonism is a phenomenon in which binding of different agonists to a target GPCR promotes distinct receptor conformations that bias cellular signaling toward and away from a subset of pathways, e.g., activation of different G protein or arrestin [100][101][102][103]. Biased agonism was first introduced by Jarpe et al. [103] and has been targeted for developing selective therapeutics of GPCRs [102,104]. The first A1AR biased agonist LUF5589 was identified by assessing A1AR bias with 800 A1AR agonists and antagonists AR-G s complex systems regarding the agonist RMSD relative to the cryo-EM conformation and AR:NPxxY-G:α5 distance. The white triangles indicate the cryo-EM or simulation starting structures. Summary of specific AR-G protein interactions: (E) the ADO-bound A 1 AR prefers to bind the G i protein to the G s . The latter could not stabilize agonist ADO binding in the A 1 AR and tended to dissociate from the receptor. (F) The A 2A AR could bind both the G s and G i proteins, which adopted distinct conformations in the complexes. Adapted from reference [99] with permission from American Chemical Society. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Allosteric Modulation of Adenosine Receptors
Due to high similarity of the orthosteric site among different ARs, agonists have failed clinical trials due to off-target side effects. To overcome this problem, allosteric modulators have been developed that bind topographically distinct "allosteric" sites in the receptor ( Figure 3) and regulate the effects of orthosteric ligands. PAMs potentiate the effects of orthosteric ligands, whereas negative allosteric modulators (NAMs) do the opposite. AR allosteric modulators provide subtype selectivity reducing the off-target side effects and hence represent a promising approach to selective drug design. Gao et al. [111] characterized the binding and functional antagonism of fluorescent conjugates of xanthine amine congener (XAC) and SCH442416 to the A2AAR. Among antagonists tested, MRS7322, MRS7396, MRS7416, XAC245 and XAC630 behaved as allosteric antagonists of A2AAR, whilst MRS7395 and XAC488 acted as competitive antagonists [111]. Allosteric antagonists MRS7396 and 5-(N,N-hexamethylene)amiloride (HMA) were more potent than MRS7416 in displacing [ 3 H]ZM241385 binding, while MRS7396, XAC630 and HMA were less potent than radioligand binding in displacing MRS7416 binding [111]. Mutation

Allosteric Modulation of Adenosine Receptors
Due to high similarity of the orthosteric site among different ARs, agonists have failed clinical trials due to off-target side effects. To overcome this problem, allosteric modulators have been developed that bind topographically distinct "allosteric" sites in the receptor ( Figure 3) and regulate the effects of orthosteric ligands. PAMs potentiate the effects of orthosteric ligands, whereas negative allosteric modulators (NAMs) do the opposite. AR allosteric modulators provide subtype selectivity reducing the off-target side effects and hence represent a promising approach to selective drug design. Gao et al. [111] characterized the binding and functional antagonism of fluorescent conjugates of xanthine amine congener (XAC) and SCH442416 to the A 2A AR. Among antagonists tested, MRS7322, MRS7396, MRS7416, XAC245 and XAC630 behaved as allosteric antagonists of A 2A AR, whilst MRS7395 and XAC488 acted as competitive antagonists [111]. Allosteric antagonists MRS7396 and 5-(N,N-hexamethylene)amiloride (HMA) were more potent than MRS7416 in displacing [ 3 H]ZM241385 binding, while MRS7396, XAC630 and HMA were less potent than radioligand binding in displacing MRS7416 binding [111]. Mutation of D52N in the sodium site of A 2A AR changed the affinity of HMA and MRS7396, indicating preferences for different A 2A AR conformations [111].
In one study, we performed GaMD simulations to identify binding modes of two prototypical PAMs in the A 1 AR [112]. VCP171 and PD81723 were initially placed around 20 Å away from the receptor. Amber [113] and NAMD [114] software packages were used to perform enhanced sampling of PAMs binding to the receptor. GaMD predicted the binding modes of both the PAMs near ECL2 site which was consistent with the experimental mutagenesis results [115]. In the PAM-bound state of the receptor, the NECA agonist exhibited lower fluctuations in the orthosteric site. In contrast, without the PAM, the agonist was observed to be very dynamic in the orthosteric site and could even dissociate in some of the simulations.
Very recently, using the first cryo-EM structure of A 1 AR bound to a PAM MIPS521 in presence of bound ADO agonist and G protein, we performed further GaMD simulations to characterize molecular basis of allosteric modulation [39]. GaMD simulations showed that the MIPS521 PAM molecule could stabilize the ADO in the orthosteric pocket ( Figure 4). Even in the absence of G protein, the PAM molecule could show positive cooperativity stabilizing the ADO agonist. We could observe deactivation during GaMD simulation of A 1 AR without the PAM and G protein, whereas the presence of MIP521 could slow deactivation of A 1 AR without the G protein bound (Figure 4). This showed that the MIP521 PAM could stabilize the ADO-bound of A 1 AR in a "G-protein-bound-like" conformation. of D52N in the sodium site of A2AAR changed the affinity of HMA and MRS7396, indicating preferences for different A2AAR conformations [111].
In one study, we performed GaMD simulations to identify binding modes of two prototypical PAMs in the A1AR [112]. VCP171 and PD81723 were initially placed around 20 Å away from the receptor. Amber [113] and NAMD [114] software packages were used to perform enhanced sampling of PAMs binding to the receptor. GaMD predicted the binding modes of both the PAMs near ECL2 site which was consistent with the experimental mutagenesis results [115]. In the PAM-bound state of the receptor, the NECA agonist exhibited lower fluctuations in the orthosteric site. In contrast, without the PAM, the agonist was observed to be very dynamic in the orthosteric site and could even dissociate in some of the simulations.
Very recently, using the first cryo-EM structure of A1AR bound to a PAM MIPS521 in presence of bound ADO agonist and G protein, we performed further GaMD simulations to characterize molecular basis of allosteric modulation [39]. GaMD simulations showed that the MIPS521 PAM molecule could stabilize the ADO in the orthosteric pocket ( Figure 4). Even in the absence of G protein, the PAM molecule could show positive cooperativity stabilizing the ADO agonist. We could observe deactivation during GaMD simulation of A1AR without the PAM and G protein, whereas the presence of MIP521 could slow deactivation of A1AR without the G protein bound (Figure 4). This showed that the MIP521 PAM could stabilize the ADO-bound of A1AR in a "G-protein-boundlike" conformation. (e-h) Distance between the intracellular ends of TM3 and TM6 (measured as the distance in Å between charge centers of residues R 3.50 and E 6.30 ) in the absence (e) or presence (f) of MIPS521, Gi2 (g) or both (h). Each condition represents three GaMD simulations, with each simulation trace displayed in a different color (black, red, blue). The lines depict the running average over 2 ns. Reprinted from reference [39] with permission from Springer Nature.
SuMD simulations were also performed to investigate the binding of LUF6000 PAM to the A3AR [116]. The SuMD simulations characterized the binding pathway of LUF6000 PAM to the A3AR in presence of agonist ADO. In presence of the agonist, LUF6000 bound to the receptor in two different poses. First, the PAM bound to the ECL2 site changing its conformation which further induced energetically favorable agonist interactions in the (e-h) Distance between the intracellular ends of TM3 and TM6 (measured as the distance in Å between charge centers of residues R 3.50 and E 6.30 ) in the absence (e) or presence (f) of MIPS521, G i2 (g) or both (h). Each condition represents three GaMD simulations, with each simulation trace displayed in a different color (black, red, blue). The lines depict the running average over 2 ns. Reprinted from reference [39] with permission from Springer Nature.
SuMD simulations were also performed to investigate the binding of LUF6000 PAM to the A 3 AR [116]. The SuMD simulations characterized the binding pathway of LUF6000 PAM to the A 3 AR in presence of agonist ADO. In presence of the agonist, LUF6000 bound to the receptor in two different poses. First, the PAM bound to the ECL2 site changing its conformation which further induced energetically favorable agonist interactions in the orthosteric site. Second, the PAM stably bound near the mouth of the orthosteric site stabilizing the ternary complex with agonist-bound receptor. In another study, SuMD simulations were performed to study the allosteric role of sodium ion binding to the A 2A AR [117]. The sodium ion has been proposed as a NAM as it was observed in a distal site compared to the orthosteric site [71]. The sodium ion was able to coordinate the inactive A 2A AR without any conformational changes in the SuMD simulations. However, the TM helices rearranged to accommodate the sodium ion in an intermediate-active state of A 2A AR. Deganutti et al. performed SuMD simulations to capture binding of the antagonist 13B, bitopic agonist VCP746 and PAMs PD81723 and VCP171 to the A 1 AR [110]. The SuMD simulations showed that PAMs can bind to several receptor sites rather than a single allosteric pocket. In absence of the NECA agonist, PAM showed partial agonism behavior. Despite having structural similarities between PAMs, different binding paths were observed in the simulations, which revealed dramatic effects of subtle chemical modifications in the ligand structure.

Ligand Binding to Adenosine Receptors
Ligand binding kinetics, especially the dissociated rate, have recently been recognized to be potentially more relevant for drug design. MD simulations have been performed in order to understand ligand binding/unbinding of ARs. In a study combining MD simulations and kinetic radioligand binding experiments, Guo et al. investigated the dissociation mechanism of an antagonist ZM241385 to the A 2A AR [118]. MD simulations that captured the dissociation of the antagonist ZM241385 from the A 2A AR helped identify the residues in the ligand unbinding pathway. Experiments validated that mutation of these residues could influence ligand's dissociation rate dramatically even though the binding affinity was barely changed. This study demonstrates that receptor structural elements that are not important to binding affinity can prove key to ligand kinetics.
SuMD simulations were performed to study the binding interactions of ADO and its metabolite inosine ligand to the A 2 AR in both the active-intermediate and its G proteinbound conformations [119]. During the SuMD simulations of ADO-bound A 2A AR, the ligand was stabilized by two hydrogen bonds formed with residues N 6.55 and E169 ECL2 in the orthosteric pocket. Conversely, inosine could form only one hydrogen bond with N 6.55 . Interestingly, ligand binding in the orthosteric pocket of both systems was remarkably influenced by the presence of G protein. In another SuMD study, Sabbadin et al. explored the recognition pathway of ADO by the A 2A AR [120]. SuMD simulations showed that ECL3 represents a possible metastable binding site for ADO. During the binding process, ECL3 helped orient the adenosine ribose ring toward orthosteric entrance. In the orthosteric binding site, ribose moiety of the ligand experienced dynamic flipping between "ribosedown" and "ribose-up" conformations. Bolcato et al. performed SuMD simulations to study binding of subtype selective antagonists LC4 and Z48 binding to A 1 AR and A 2A AR, respectively [121]. The simulations showed that receptor-ligand recognitions were multistep processes involving intermediate states to guide the (un)binding events. Overall, Z48 favored A 2A AR over A 1 AR, forming classic antagonist fingerprint interactions at the orthosteric site. In case of the A 1 AR, a water molecule was seen playing a key role in stabilizing LC4 at the orthosteric pocket, which was not observed in the case of A 2A AR. In another study, Deganutti et al. [122] performed SuMD simulations to investigate the binding/unbinding pathways of five different A 1 AR agonists including ADO, CPA, NECA, HOCPA and BnOCPA. The SuMD simulations showed that the ligand followed the binding paths involving mainly ECL2, the top part of TM1, TM2, TM6 and TM7. The ligand dissociated following similar paths; however, ECL2 was less engaged. These pathways were further supported by alanine mutagenesis experiments.
Recently, we performed GaMD simulations to determine the binding and dissociation pathways of caffeine (CFF) antagonist to the A 2 AR [123]. The X-ray structure of A 2A AR in complex with CFF (PDB: 5MZP) [36] was used as initial conformation. A total of 10 CFF molecules were placed randomly at a distance >15 Å from the extracellular surface of the A 2A AR. Spontaneous binding and dissociation of CFF in the receptor were successfully captured in the GaMD simulations. A main binding pathway of CFF to the A 2A AR was identified from the 63-ns GaMD equilibration ( Figure 5A). CFF reached its binding site of the A 2A AR through interacting with ECL2, ECL3, TM7 and finally the receptor orthosteric site ( Figure 5D). GaMD production simulations captured a slightly different binding pathway when the orthosteric pocket was already occupied by one CFF molecule ( Figure 5B). The second CFF explored a region between ECL3 and TM7 during this binding process ( Figure 5E). The dissociation pathway of CFF was mostly the reverse of the dominant binding pathway (Figure 5C,F).
Molecules 2022, 27, x FOR PEER REVIEW 13 of 27 of 10 CFF molecules were placed randomly at a distance >15 Å from the extracellular surface of the A2AAR. Spontaneous binding and dissociation of CFF in the receptor were successfully captured in the GaMD simulations. A main binding pathway of CFF to the A2AAR was identified from the 63-ns GaMD equilibration ( Figure 5A). CFF reached its binding site of the A2AAR through interacting with ECL2, ECL3, TM7 and finally the receptor orthosteric site ( Figure 5D). GaMD production simulations captured a slightly different binding pathway when the orthosteric pocket was already occupied by one CFF molecule ( Figure 5B). The second CFF explored a region between ECL3 and TM7 during this binding process ( Figure 5E). The dissociation pathway of CFF was mostly the reverse of the dominant binding pathway ( Figure 5C,F).

Lipid Interactions with Adenosine Receptors
Lipid bilayers have been shown to modulate GPCR functions, including conformation stability, ligand binding and oligomerization [124][125][126]. Modulatory effects mediated via changes in the physical properties of membrane, such as thickness, curvature and surface tension, have been extensively studied using experimental and computational methods [126,127]. Here, we discuss MD studies focusing on AR-lipid interactions.
In one study, GaMD simulations were applied to study the relationship between the lipid environment and A1AR activation state [128]. The cryo-EM structure of the active ADO-bound A1AR coupled with the Gi protein (A1AR-Gi) [38] and the X-ray structure [35,36] of the inactive antagonist PSB36-bound A1AR (PSB36-A1AR) were used to perform GaMD simulations. The A1AR was embedded in a 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine (POPC) lipid bilayer. GaMD simulations suggested that the membrane lipids play a critical role in stabilizing different states of the A1AR. Different structural flexibility profiles were identified in GaMD simulations of the inactive and active A1AR. Compared with the inactive A1AR, higher fluctuations were found at the ECL2 region, intracellular ends of TM5 and TM6 in the active A1AR. The receptor TM domain and the Reprinted from reference [123] with permission from Frontiers.

Lipid Interactions with Adenosine Receptors
Lipid bilayers have been shown to modulate GPCR functions, including conformation stability, ligand binding and oligomerization [124][125][126]. Modulatory effects mediated via changes in the physical properties of membrane, such as thickness, curvature and surface tension, have been extensively studied using experimental and computational methods [126,127]. Here, we discuss MD studies focusing on AR-lipid interactions.
In one study, GaMD simulations were applied to study the relationship between the lipid environment and A 1 AR activation state [128]. The cryo-EM structure of the active ADO-bound A 1 AR coupled with the G i protein (A 1 AR-Gi) [38] and the X-ray structure [35,36] of the inactive antagonist PSB36-bound A 1 AR (PSB36-A 1 AR) were used to perform GaMD simulations. The A 1 AR was embedded in a 1-palmitoyl-2-oleoyl-glycero-3phosphocholine (POPC) lipid bilayer. GaMD simulations suggested that the membrane lipids play a critical role in stabilizing different states of the A 1 AR. Different structural flexibility profiles were identified in GaMD simulations of the inactive and active A 1 AR. Compared with the inactive A 1 AR, higher fluctuations were found at the ECL2 region, intracellular ends of TM5 and TM6 in the active A 1 AR. The receptor TM domain and the ligands were rather rigid. However, the G protein coupled to the active A 1 AR exhibited high flexibility, especially in the, α4-β5 loop and α4-β6 loop and α5 helix of the Gα subunit and terminal ends of the G βγ subunits.
The -S CD order parameter values obtained from GaMD simulations agreed well with experimental data. In NMR experiments, the -S CD for the fifth carbon C-H bond of POPC was at~0.18-0.20 [129]. The -S CD value of POPC's fifth carbon atom was~0.17 ± 0.02 in the lower leaflet in the active A 1 AR system. It increased to~0.20 ± 0.02 in the inactive A 1 AR system ( Figure 6). The -S CD value of the ninth carbon C-H bond in POPC calculated from GaMD simulations was~0.10, being consistent with the NMR experiments [129]. Additionally, the GaMD simulations suggested that POPC lipids in the lower leaflet of the inactive A 1 AR system were less fluid than in the active A 1 AR system. The similar -S CD values of sn-2 acyl chains of POPC molecules in the upper leaflet were identified in the inactive and active A 1 AR systems. However, the -S CD for the lower leaflet in the inactive A 1 AR system was larger than those in the active A 1 AR system. This finding correlated well with the TM6 outward movement in the active A 1 AR, which caused higher inclination of the C-H bonds to be aligned along the bilayer normal. In summary, GaMD simulations have revealed strongly coupled dynamics between a GPCR and the membrane lipids that depend on the receptor activation state. ligands were rather rigid. However, the G protein coupled to the active A1AR exhibited high flexibility, especially in the, α4-β5 loop and α4-β6 loop and α5 helix of the Gα subunit and terminal ends of the Gβγ subunits. The -SCD order parameter values obtained from GaMD simulations agreed well with experimental data. In NMR experiments, the -SCD for the fifth carbon C-H bond of POPC was at ~0.18-0.20 [129]. The -SCD value of POPC's fifth carbon atom was ~0.17 ± 0.02 in the lower leaflet in the active A1AR system. It increased to ~0.20 ± 0.02 in the inactive A1AR system (Figure 6). The -SCD value of the ninth carbon C-H bond in POPC calculated from GaMD simulations was ~0.10, being consistent with the NMR experiments [129]. Additionally, the GaMD simulations suggested that POPC lipids in the lower leaflet of the inactive A1AR system were less fluid than in the active A1AR system. The similar -SCD values of sn-2 acyl chains of POPC molecules in the upper leaflet were identified in the inactive and active A1AR systems. However, the -SCD for the lower leaflet in the inactive A1AR system was larger than those in the active A1AR system. This finding correlated well with the TM6 outward movement in the active A1AR, which caused higher inclination of the C-H bonds to be aligned along the bilayer normal. In summary, GaMD simulations have revealed strongly coupled dynamics between a GPCR and the membrane lipids that depend on the receptor activation state. Figure 6. Differentiated order parameters of lipid molecules were found in simulation systems of the inactive and active A1AR: (A) inactive A1AR using dihedral-boost GaMD, (B) active A1AR using dihedral-boost GaMD, (C) inactive A1AR using dual-boost GaMD and (D) active A1AR using dualboost GaMD. Red diamond lines represent the average -SCD order parameters for the cytoplasmic lower leaflet and blue diamond lines for the extracellular upper leaflet. Reprinted from reference [128] with permission from Wiley. Figure 6. Differentiated order parameters of lipid molecules were found in simulation systems of the inactive and active A 1 AR: (A) inactive A 1 AR using dihedral-boost GaMD, (B) active A 1 AR using dihedral-boost GaMD, (C) inactive A 1 AR using dual-boost GaMD and (D) active A 1 AR using dualboost GaMD. Red diamond lines represent the average -S CD order parameters for the cytoplasmic lower leaflet and blue diamond lines for the extracellular upper leaflet. Reprinted from reference [128] with permission from Wiley. MD simulations were also performed to explore the effects of different lipids in antagonist caffeine binding to A 2A AR [130]. Only POPC, a mixture of POPC and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylethanolamine (POPE), and cholesterol rich lipid environment were used in the study. Simulations showed that H8 folding depended on the lipid environment. Cholesterol was observed to bind a cleft between TM1 and TM2, which stabilized a distinct caffeine binding pose. This further highlighted the importance of using physiological cholesterol concentration in MD simulations. Recently, Bruzzese et al. performed MD simulations to study the role of different lipids in activation of A 2A AR [131]. Inactive A 2A AR in the apo or agonist (ADO or NECA) bound conformations were used as starting structures in either 1,2-dioleoyl-sn-glycerol-3-phosphoglycerol (DOPG) or 1,2-dioleoyl-sn-glycerol-3-phosphocholine (DOPC) lipid bilayer. DOPC could facilitate the transition of the inactive A 2A AR to an intermediate conformation. In DOPG with apo conformation of A 2A AR, the receptor could also transit to an intermediate conformation.
Interestingly, the agonist-bound A 2A AR transitioned into a fully active state. These variations were attributed to the allosteric effects mediated by the lipid and presence/absence of an agonist. Leonard et al. performed 35 µs MD simulations to study the preference for lipid solvation in the A 2A AR [132]. Free energy of lipid solvation was calculated for different activation states of A 2A AR with different lipids. The results showed that the inactive state preferred unsaturated lipids over the saturated ones for formation of the first solvation lipid shell around the receptor. This preference enhanced even more in the partially active A 2A AR state as compared to that of the inactive receptor.
In addition to the above-mentioned all-atom MD simulations, coarse-grained models were also widely used to reach longer time scale and/or simulate larger systems. For example, both all-atom and coarse-grained MD simulations were performed to study the effects of membrane lipids and cations in the inactive, intermediate and fully active G protein-bound conformations of A 2A AR [133]. The study was performed in both detergent micelles and lipid bilayer environment. MD simulations suggested that phosphatidylinositol bisphosphate 2 (PIP2) interacted with the A 2A AR intracellular residues, which could reduce the flexibility of the receptor in the inactive state and limit the transition to the active-intermediate state. For the fully active A 2A AR, PIP2 stabilized the receptor-G protein complex. However, such stabilizing interactions were absent in the non-ionic micelles. The level of activation microswitches observed in the POPC lipid bilayer and detergent micelles were different, suggesting a rheostat model of GPCR activation microswitches as compared to a binary switch model. Song et al. [134] used coarse-grained MD simulations to study A 2A AR-lipid interactions. Simulations showed that different kinds of lipids could interact with nine binding sites on the receptor. The lipids were observed to allosterically modulate activation of A 2A AR. PIP2 could help stabilize the active state of the A 2A AR by stabilizing outward movement of TM6 and enhancing the AR-G protein interactions. These studies strengthened the notion that lipids affect GPCR signaling and function by allosterically modulating its membrane lipid environment, being consistent with experimental findings [135]. For example, Huang et al. [135] identified cholesterol as a weak PAM of the A 2A AR by combining 9 F NMR, computational analysis and G protein assays. Their GTP hydrolysis assays showed a marginal increase in basal activity with increasing cholesterol, in addition to a weak enhancement in the agonist potency. Furthermore, their 19 F NMR data suggested that the enhancement resulted from an increase in the receptor's active state population and a G protein-bound pre-coupled state.

Computer-Aided Drug Discovery of Adenosine Receptors
More than 515 million compounds are available in the ZINC chemical database [136]. It is impossible to perform high-throughput screening of all the chemical compounds. In this regard, virtual screening is a valuable strategy to efficiently screen a large number of compounds and select only a small subset of promising compounds for testing in in vitro and in vivo experiments [31,137,138]. Rodríguez et al. [139] identified 20 potent ligands by molecular docking from a pool of more than 6.7 million commercially available com-pounds [139]. MD simulations have been incorporated into virtual screening to increase its success rate by providing a more reliable receptor structural ensemble and more accurate binding free energy calculation. In the following, we discuss novel binding site identification and drug design of ARs.

Drug Binding Sites of Adenosine Receptors
The orthosteric site has long been targeted to design drugs for ARs. However, this site is usually conserved across four subtypes of ARs. Consequently, it is critical to identify novel, less conserved binding sites to design selective ligands of ARs. Caliman et al. [140] used FTPMap to identify potential binding sites of the A 2A AR in the 20 A 2A AR crystal structures and 30 receptor clusters with 1.57 µs and 1.75 µs MD simulations of the apo receptor generated from structures of 3EML and 3QAK. Five "non-orthosteric" sites were identified, including the extracellular region of TM3 and TM4, lipid interface of TM5 and TM6, intracellular region between TM1 and TM7, G protein binding site among TM2, TM3, TM6 and TM7 and intracellular region of TM3 and TM4 [140].
Another method, the site-identification by ligand competitive saturation (SILCS) [141,142] method, was also developed to predict the location and approximate affinities of small molecular fragments on a target receptor surface by performing MD simulations. Moreover, MD simulations and mutation experiments were combined to identify residue hotspots of ARs for ligand binding. Wang et al. [143] explored the importance of specific residue hotspots by mutating Leu 6.51 to Val specific to the A 2B AR in the A 2A AR. They incorporated MD simulations and radioligand binding assays to validate the mutation. A selective A 2A AR antagonist ZM241385 indicated decreased affinity for the L249V 6.51 mutant of A 2A AR, suggesting the important role of this Leucine residue.

Binding Free Energy Calculations of Adenosine Receptors
One of the fundamental assumptions of ligand-based drug design is that similar molecules will have similar biological activity. Such an approach often fails when a small change in the ligand structure leads to a drastic difference in binding affinity [144]. In this regard, MD simulation and binding free energy calculation methods, including molecular mechanics Poisson-Boltzmann surface area (MM/PBSA), molecular mechanics generalized Born surface area (MM/GBSA) and free energy perturbation (FEP), have become valuable tools to predict ligand binding affinities for drug design [145,146]. MM/PBSA and MM/GBSA are popular methods for binding free energy prediction in drug design since they are more accurate than most scoring functions of molecular docking [146]. For example, Lenselink et al. [147] identified two ligands of A 2A AR by applying MM/GBSA to rescore binding poses from Glide docking. Compared with MM/PBSA or MM/GBSA, MD/FEP is a more accurate binding free energy calculation method. Therefore, MD/FEP has been routinely used in drug design projects [45,145,[148][149][150][151][152]. For the A 2A AR, MD/FEP calculations suggested the loss of binding of 3-deazaadenosine due to modification of a single heavy atom in ADO [150]. Even for such a small modification, MD/FEP was able to successfully predict the change of ligand binding affinity. One lead compound predicted from MD/FEP was synthesized and experimentally verified to be a full agonist equipotent to ADO. Furthermore, MD/FEP suggested that water in the binding site provided a major driving force for the ligand-protein interaction.
Fragment-based drug discovery relies on successful optimization of weak ligands for binding affinity and selectivity. MD/FEP has been widely applied in the fragment-to-lead optimization process. Matricon et al. [149] discovered that the benzothiazole fragment bound to the A 1 AR through hydrogen bond interactions with N 6.55 . Through a single vector growth, they designed nine additional compounds from the initial fragment, all of which showed improved binding affinity as suggested by MD/FEP calculations [149]. Eventually, the compounds were synthesized and their binding affinities in the A 1 AR were experimentally measured by radioligand binding assays [149]. The resulting ligands led to >1000-fold improvements of binding affinity and nearly 40-fold higher subtype selectivity. Furthermore, MD/FEP has been used to identify the ligand binding pose. In a proof-ofconcept study, Jespers et al. [45] presented a robust protocol based on iterations of MD/FEP calculations, chemical synthesis, biophysical mapping and structural biology to determine the binding pose of a series of antagonists in the A 2A AR. Eight binding site mutations in the A 2A AR were initially analyzed with MD/FEP calculations. The suggested ligand binding pose was subsequently used to guide the design of the new analogues. Remarkably, the binding affinities obtained from MD/FEP prediction were highly consistent with the experimental data. Furthermore, the predicted binding poses were highly consistent with the experimental structures.

Design of Biased Agonists of Adenosine Receptors
Wall et al. [106] discovered biased agonist BnOCPA of A 1 AR that favors G ob over other G protein subtypes and hence produces analgesic effects without sedation, bradycardia and hypotension. MD simulations revealed the BnOCPA binding modes in the receptor. BnOCPA binding resulted in unique active-and inactive-like conformations of the receptor during the simulations. Similarly, Baltos et al. [153] investigated the A 3 AR biased agonists and discovered that several methyluronamide nucleoside derivatives exhibited biased agonism. In particular, molecular docking of MRS5679, an (N)-methanocarba nucleoside derivative with an extended C2 substituent, to a homology model of the A 3 AR provided important insights into the molecular mechanism of biased signaling. Binding of the MRS5679 biased agonist favored a distinct conformation of the A 3 AR with outward movement of the TM2 extracellular domain. Valant et al. [109] combined the endogenous agonist and allosteric modulator to design a bitopic ligand VCP746 as a biased agonist of A 1 AR that favored cAMP pathway over ERK1/2 phosphorylation pathway. Molecular docking calculations and SuMD simulations [110,154] were used to study ligand binding to the receptor (see Section 3.4).

Design of Allosteric Modulators of Adenosine Receptors
Virtual screening has been widely used for agonist/antagonist design targeting GPCRs [31]. However, it is rather challenging to apply virtual screening to discover allosteric modulators with high potency. This is largely due to the lower affinity of allosteric modulators compared with the agonists/antagonists and high receptor flexibility of the target sites. Induced-fit docking [155] and ensemble docking [156] have been developed to account for receptor flexibility in virtual screening. The structural ensembles of target receptors could be taken from NMR structures, homology modeling, MD simulations or any other conformational sampling method [156][157][158][159]. It has been applied in AR allosteric modulator design. The receptor structural ensemble is often generated from MD simulations. Particularly, enhanced sampling MD simulations could sample larger conformational space of the drug target and have thus been shown to increase the success rate of finding allosteric modulators of the M 2 muscarinic GPCR [160]. In a recent study [161], we tested whether the receptor ensemble generated from GaMD simulations could increase the docking accuracy of discovering PAMs in the A 1 AR. Extensive retrospective ensemble docking calculations of PAM binding to the A 1 AR were performed using GaMD simulations and AutoDock (Figure 7) [162]. The dihedral and dual boost GaMD implemented in the AMBER and NAMD were performed. The rigid-body docking at short, medium and long levels and flexible docking were all evaluated in the ensemble docking protocol. The boost potential obtained for the same system with AMBER was larger than with NAMD in the GaMD simulations, which is due to different algorithms used to calculate the system potential statistics [112]. Accordingly, larger conformation space of the receptor was sampled in the AMBER simulations, leading to improved docking performance. Correction of docking score by the GaMD reweighted free energy of the receptor structural cluster further improved the docking performance. With GaMD reweighted scores, ranking by the average binding energy (BE avg ) performed better than by the minimum binding energy (BE min ) in terms of the area under the receiver operating characteristic curves (AUC), enrichment factors (EFs) metrics. The receptor ensemble obtained from AMBER dual boost GaMD simulations of the VCP171 PAM-bound ADO-A 1 AR-Gi outperformed other receptor ensembles for docking. This ensemble consisted of conformations of the holo A 1 AR with PAM bound at the ECL2 allosteric site. Interactions between the PAM and receptor ECL2 might induce more reliable conformations for PAM binding, which were otherwise difficult to sample in the simulations of PAM-free (apo) A 1 AR. Ensembles obtained from dual boost GaMD performed better than the dihedral-boost GaMD for docking, suggesting that GaMD with a higher acceleration level was needed to sufficiently sample conformational space of the GPCR PAM binding site. docking score by the GaMD reweighted free energy of the receptor structural cluster further improved the docking performance. With GaMD reweighted scores, ranking by the average binding energy ( ) performed better than by the minimum binding energy ( ) in terms of the area under the receiver operating characteristic curves (AUC), enrichment factors (EFs) metrics. The receptor ensemble obtained from AMBER dual boost GaMD simulations of the VCP171 PAM-bound ADO-A1AR-Gi outperformed other receptor ensembles for docking. This ensemble consisted of conformations of the holo A1AR with PAM bound at the ECL2 allosteric site. Interactions between the PAM and receptor ECL2 might induce more reliable conformations for PAM binding, which were otherwise difficult to sample in the simulations of PAM-free (apo) A1AR. Ensembles obtained from dual boost GaMD performed better than the dihedral-boost GaMD for docking, suggesting that GaMD with a higher acceleration level was needed to sufficiently sample conformational space of the GPCR PAM binding site. , GaMD simulations were carried out to construct structural ensembles to account for the receptor flexibility. Meanwhile, a compound library was prepared for 25 known PAMs of the A1AR and 2475 decoys obtained from the DUD-E with openbabel 2.4.1. Ensemble docking was then performed to identify the PAMs for which the AUC and enrichment factors were calculated to evaluate docking performance. Both rigid-body and flexible docking were tested using AutoDock. Reprinted from reference [161] with permission from Elsevier.
Overall, flexible docking performed significantly better than the rigid-body docking at different levels with AutoDock. This suggested that the flexibility of protein side chains in ensemble docking also played an important role. The side chains of representative re- in the A 1 AR: Starting from the cryo-EM structure of the active ADO-Gi-bound A 1 AR (6D9H) and docking model of PAM VCP171-bound A 1 AR (ADO-A 1 AR-Gi-VCP171), GaMD simulations were carried out to construct structural ensembles to account for the receptor flexibility. Meanwhile, a compound library was prepared for 25 known PAMs of the A 1 AR and 2475 decoys obtained from the DUD-E with openbabel 2.4.1. Ensemble docking was then performed to identify the PAMs for which the AUC and enrichment factors were calculated to evaluate docking performance. Both rigid-body and flexible docking were tested using AutoDock. Reprinted from reference [161] with permission from Elsevier.
Overall, flexible docking performed significantly better than the rigid-body docking at different levels with AutoDock. This suggested that the flexibility of protein side chains in ensemble docking also played an important role. The side chains of representative receptor structures obtained from GaMD simulations might be still in unfavored conformations for PAM binding. Flexible docking of protein side chains could then alleviate this problem to achieve better performance. In summary, GaMD simulations and flexible docking greatly improved the docking performance by effectively accounting for the protein flexibility in both the backbone and side chains. Such an ensemble docking protocol will greatly facilitate future allosteric drug design of the A 1 AR.

Design of Bitopic Ligands of Adenosine Receptors
Bitopic ligands contain hybrid molecular structures with pharmacophores of both orthosteric and allosteric ligands. Bitopic drug design has gained popularity because of its multivariate functions including increased efficacy and biased agonism in one ligand. Narlawar et al. [163] designed bitopic ligands of A 1 AR combining orthosteric and allosteric pharmacophores with increasing length of linker between them. In particular, LUF6258 with a 9-carbon atom spacer between the allosteric and orthosteric moieties stood out with increased efficacy. The homology model of the A 1 AR was used for docking of compounds and the binding poses were analyzed. LUF6258 had its adenosine part in the orthosteric pocket, whereas the allosteric moiety could extend out to the extracellular space and interacted with the ECL2 region. Valant et al. [109] rationally designed a new bitopic ligand VCP746 by combining endogenous agonist adenosine and PAM VCP171. VCP746, in an A 1 AR-mediated inhibition assay, showed 30-fold biased agonism toward cAMP pathway in comparison to ERK1/2 phosphorylation pathway. In a follow up study [154], structural changes were made to VCP746, and structural-activity relationship experiments suggested that the 4-(trifluoromethyl) phenyl allosteric pharmacophore group was responsible for the ligand bias effect in A 1 AR. Docking calculations using the A 1 AR homology model suggested that the allosteric and linker region of the bitopic ligand explored the extracellular vestibule of the receptor. Deganutti et al. [110] applied SuMD to simulate the binding of VCP746 to A 1 AR and proposed the binding modes of the bitopic ligand. The adenosine moiety of the ligand bound stably at the orthosteric site, whereas the VCP171 part bound near the ECL2 region. These molecular details further helped understanding the functional mechanism of the bitopic ligand. These studies provide evidence that combining orthosteric and allosteric pharmacophores improves pharmacological drug properties such as on-target efficacy, biased agonism, less on-target side effects, etc.

Conclusions
ARs have served as established drug targets. Computer-aided structure-based drug design approach has proven useful as while traditional agonists and antagonists often induce adverse side effects, new paradigms have emerged to search for novel biased agonists and allosteric modulators that can serve as selective drug molecules of the ARs. In this regard, because the target sites of these new ligands often involve regions on the receptor surface, it is crucial to account for receptor flexibility in order to design potent drug molecules. Additionally, kinetics of ligand binding has recently been recognized to be potentially more relevant for drug design. In particular, the dissociation rate constant that determines the drug residence time appears to better correlate with drug efficacy than the binding free energy [164,165]. However, ligand kinetic rates are even more difficult to compute than the binding free energies, largely due to slow processes of ligand binding and dissociation over long time scales [165]. Remarkable advances in MD approaches and computer hardware have paved the way to capture the ligand binding/unbinding processes of ARs in molecular details [166], which is expected to play a more important role in drug design in the future. In this context, advanced computational platforms (e.g., ANTON 3 [167] and Deep Docking [168]) and improved simulation techniques, including the coarse-grained MD, FEP, ligand and peptide GaMD [169,170], SuMD [110], Metadynamics [93], machine learning [171] and deep learning [168,172,173], will continue to drive rational computer-aided design of more potent and selective drug molecules of ARs.

Conflicts of Interest:
The authors declare no conflict of interest.